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In this paper, we investigate the emergence of a predator-prey model with Beddington-DeAngelis- 
type functional response and reaction-diffusion. We derive the conditions for Hopf and Turing 
bifurcation on the spatial domain. Based on the stability and bifurcation analysis, we give the 
, spatial pattern formation via numerical simulation, i.e., the evolution process of the model near the 

coexistence equilibrium point. We find that for the model we consider, pure Turing instability gives 
1^-^ ' birth to the spotted pattern, pure Hopf instability gives birth to the spiral wave pattern, and both 

Hopf and Turing instability give birth to stripe-like pattern. Our results show that reaction-diffusion 
i i model is an appropriate tool for investigating fundamental mechanism of complex spatiotemporal 

PQ ' dynamics. It will be useful for studying the dynamic complexity of ecosystems. 
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Si . 1. INTRODUCTION 

C3 i 

Mathematical models have played an important role throughout the history of ecology. Models in ecology serve 
a variety of purposes, which range from illustrating an idea to parameterizing a complex real-world situation. They 
are used to make general predictions, to guide management practices, and to provide a basis for the development of 
statistical tools and testable hypotheses (ToL 4fl 51 1 . 
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A fundamental goal of theoretical ecology is to understand how the interactions of individual organisms with each 
other and with the environment determine the distribution of populations and the structure of communities 0] . As 
we know, our ecological environment is a huge and highly complex system. This complexity arises in part from 



the diversity of biological species, and also from the complexity of every individual organism 2^, 3^]. The dynamic 
behavior of predator-prey model has long been and will continue to be one of the dominant themes in both ecology 
and mathematical ecology due to its universal existence and importance 0, • 

Before the 1970s, ecological population models typically used ordinary differential equations, seeking equilibria and 
analyzing stability. The early models provided important insights, such as when species can stably coexist and when 
predator and prey density oscillate over time [loj ]. 
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We live in a spatial world, and spatial patterns are ubiquitous in nature, which modify the temporal dynamics and 
stability properties of population density at a range of spatial scales, whose effects must be incorporated in temporal 
ecological models that do not represent space explicitly. And the spatial component of ecological interactions has 
been identified as an important factor in how ecological communities are shaped (Io| . Empirical evidence suggests 
that the spatial scale and structure of the environment can influence population interactions and the composition of 
communities [§]. In recent decades, the role of spatial effects in maintaining biodiversity has received a great deal of 

attention in the literature on conservation d m m, h 0, m m m . 

The past investigations have revealed that spatial inhomogcneities like the inhomogeneous distribution of nutrients 
as well as interactions on spatial scales like migration can have an important impact on the dynamics of ecological 
populations [Hoi . [Hlj . In particular it has been shown that spatial inhomogeneities promote the persistence of ecological 
populations, play an important role in speciation and stabilize population levels Spatial ecology today is still 
dominated by theoretical investigations, and empirical studies that explore the role of space are becoming more 
common due to technological advances that allow the recording of exact spatial locations [l0| • 

In 1952, Alan Turing, one of the key scientists of 20th century, mathematically showed that a system of coupled 
reaction-diffusion equations could give rise to spatial concentration patterns of a fixed characteristic length from an 
arbitrary initial configuration due to diffusion-driven instability [611 ]. The work by Turing belongs to the field of 
pattern formation, a subficld of mathematical biology. Pattern formation in nonlinear complex systems is one of 
the central problems of the natural, social, and technological sciences. The occurrence of multiple steady states and 
transitions from one to another after critical fluctuations, the phenomena of excitability, oscillations, waves and the 
emergence of macroscopic order from microscopic interactions in various nonlinear noncquilibrium systems in nature 
and society have been the subject of many theoretical and experimental studies [50j]. It has been qualitatively shown 
that Turing models can indeed imitate biological patterns [37| . 

The study of biological pattern formation has gained popularity since the 1970s, and Segel and Jackson [56j were 
the first to apply Turing's ideas to a problem in population dynamics: the dissipative instability in the predator-prey 
interaction of phytoplankton and herbivorous copepods with higher herbivore motility. At the same time, Gierer and 
Meinhardt [2 11 ] gave a biologically justified formulation of a Turing model and studied its properties by employing 
numerical simulations. Levin and Segel [HI ] suggested this scenario of spatial pattern formation was a possible origin 
of planktonic patchiness. 

Spatial patterns and aggregated population distributions are common in nature and in a variety of spatiotemporal 
models with local ecological interactions 54] . And the understanding of patterns and mechanism of species spatial 
dispersal is an issue of current interest in conservation biology and ecology [3, 39, 50]. It arises from many ecological 
applications, and in p articular, p lays a major role in connection to biological invasion and epidemic spread and 
so on 0, @, [n], 0, 33, 43, 47, 43, Ho, 51 1. The field of research on pattern formation modeled by reaction-diffusion 
systems, which provides a general theoretical framework for describi ng p attern formation in systems from many diverse 
disciplines includingbiology @, 0, IMIMJisl M, M, M, M, H], chemistry [13, M, M, M, M, E3, El, M, , 



physics [12, 131|, |38|, |60| , epidemiology [43 
during the last decade. 

In general, a classical predator-prey model can be written as the form 0, 0] 



and so on, seems to be a new increasingly interesting area, particularly 



N = Nf(N) - Pg(N,P), P = h\g{N,P),P]P, 

where N and P are prey and predator densities, respectively, f(N) the prey growth rate, g(N, P) the functional 
response, e.g., the prey consumption rate by an average single predator, and h[g(N, P), P] the per capita growth rate 
of predators (also known as the "predator numerical response"), which obviously increases with the prey consumption 
rate. The most widely accepted assumption for the numerical response is the linear one Q : 

h\g(N,P),P] = eg(N,P)-ri 



where r) is a per capita predator death rate and e the conversion efficiency of food into offspring. 

In population dynamics, a functional response g(N, P) of the predator to the prey density refers to the change in 
the density of prey attached per unit time per predator as the prey density changes [jl [Hj ■ There have been several 
famous functional response types: Holling types I — III 0, [S3]; Hassell-Varley type Eqk Beddington-DeAngelis type 
by Bcddington [6j and DeAngelis et al 15] independently; the Crowley-Martin type [l3j; and the recent well-known 
ratio-dependence type by Arditi and Ginzburg Q later studied by Kuang and Beretta [35j . Of them, the Holling 
type I-III is labeled "prey-dependent" and other types that consider the interference among predators are labeled 
"predator-dependent" 0] . 
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In our previous work [64| , we studied the spatiotemporal complexity of a ratio-dependent predator-prey model with 
Michaelis-Menten-type functional response. Compare Michaelis-Menten-type functional response 

(a, h are positive constants, a capture rate and h handling time) with Bcddington-DeAngelis-type functional response 

g (N, P) = ^ (1) 

JK ' ' B + N + wP w 

(/?, B are positive constants, (3 a maximum consumption rate, B a saturation constant, w a predator interference 
parameter, a constant, w < is the case where predators benefit from cofeeding). Some scholars 0, 46|, 5t| indicat 



that where there is a small difference between the denominators, there is a world between them from the biological 
point. 

If predators do not waste time interacting with one another or if their attacks are always successful and instanta- 
neous, i.e., w — in ( [1} , then a Holling type II functional response is obtained: 

m- 0N 



B + N 



If B = in (TT|), the Michaelis-Menten-type functional response is obtained. 

Furthermore, in (59j . by comparing the statistical evidence from nineteen categories of predator- prey systems with 
three predator-dependent functional responses, Skalski and Gilliam pointed out that the predator-dependent functional 
responses could provide better descriptions of predator feeding over a range of predator-prey abundance, and in some 
cases, the Beddington-DeAngelis-type functional response performed even better [Zj|. 

Some progress has been seen in the study of predator-prey model with Beddington-DeAngelis-type functional 
response 0, |, Q [H, H [H HI, 0, [H |H IS, El . However, the research on considering reaction-diffusion to such 



a model, to our knowledge, seems rare. 

In this paper, we report a study of Turing pattern formation in a two-species reaction-diffusion predator-prey model 
with Beddington-DeAngelis-type functional response. In the next section we give a brief stability and bifurcation 
analysis of the model. Then, we present and discuss the results of numerical simulations, which is followed by the last 
section, i.e., conclusions and remarks. 



2. STABILITY AND BIFURCATION ANALYSIS 



In this paper, we mainly focus on the following predator- prey model with Beddington-DeAngelis-type functional 
response and reaction-diffusion: 

^ = r (1 - f ) N - ^§L_P + d^N, % = -r,P + d 2 ^P (2) 



where t denotes time and N, P stand for prey and predator density, respectively. All parameters are positive constants, 
r standing for maximum per capita growth rate of the prey, j3 capture rate, r\ predator death rate, w a predator 
interference parameter and K carrying capacity, which is the nonzero equilibrium population size. The diffusion 
coefficients are denoted by d\ and c?2, respectively. V 2 = + -gjp is the usual Laplacian operator in two-dimensional 
space. 

The first step in analyzing the model is to determine the equilibria (stationary states) of the non-spatial model 
obtained by setting space derivatives equal to zero, i.e., 

1 V 1 A'/ iv B+N+wP r u > B+N+wP r 'l r u - y°> 

In fact, physically, an equilibrium represents a situation without "life". It may mean no motion of a pendulum, no 
reaction in a reactor, no nerve activity, no flutter of an airfoil, no laser operation, or no circadian rhythms of biological 
clocks [57j | . And at each equilibrium point, the movement of the population dynamics vanishes. 

Eqs.® has at most three equilibria, which correspond to spatially homogeneous equilibria of the full model, Eq.([2]), 
in the positive quadrant: 

(i) (0,0) (total extinct) is a saddle point. 
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(ii) (K, 0) (extinct of the predator, or prey-only) is a stable node if e (3 < rj or if s (3 > rj and K < jgzr I a saddle if 
e [3 < T) and K > jgzr i a saddle- node if e /3 < 77 and K = jj^j ■ 

(hi) a nontrivial stationary state (N*,P*) (coexistence of prey and predator), where 

N* = — \K(rwe -s(3 + rj) + y/K 2 (rwe - e (3 + rf) 2 + 4 rKwerjB^ , 

p* _ (0e-y) jy* _ B 
wt] w ' 

To perform a linear stability analysis, we linearize the dynamic system [2] around the equilibrium point (N* , P*) for 
small space- and time-dependent fluctuations and expand them in Fourier space 

N(x, t) - N*e xt e lts , P(x, t) - P*e xt e^ s , 
and obtain the characteristic equation 

\A~k 2 D-\I\ =0, (4) 
where D = diag(di, di) and the Jacobian matrix A is given by 

9„f d p f \ ( f n fp \ 

(N*,P*) V 9n 9p J 

Now Eq.(|l]) can be solved, yielding the so called characteristic polynomial of the original problem 

A 2 -tr fe A + Afc = 0, (5) 

where 

tr fc = f n +g p - k 2 (dx + d 2 ) = tr - k 2 (d 1 + d 2 ), 

Afc = f„g p - f p g„ - k 2 (f n d 2 + g p dx) + k 4 did 2 = A - k 2 (f n d 2 + g p di) + k 4 dxd 2 . 
The roots of Eq. ([5|) yield the dispersion relation 



A = 

\d n g d p g 



A(fc) = ^(tr fe ± x /tr2-4A fc ). (6) 



We know that one type of instability (or bifurcation) will break one type of symmetry of a system, i.e., in the 
bifurcation point, two equilibrium states intersect and exchange their stability. Biologically speaking, this bifurcation 
corresponds to a smooth transition between equilibrium states. The Hopf bifurcation is space-independent and it 
breaks the temporal symmetry of a system and gives rise to oscillations that are uniform in space and periodic in 
time. The Turing bifurcation breaks spatial symmetry, leading to the formation of patterns that are stationary in 
time and oscillatory in space. 

The Hopf instability or bifurcation is an important instability in reaction-diffusion systems for which the conditions 
result in a stable limit cycle (oscillations). In terms of the linearized problem (EqHJ) the onset of Hopf instability 
corresponds to the case, when a pair of imaginary eigenvalues cross the real axis from the negative to the positive side. 
The Hopf bifurcation of a equilibrium state is reflected by a transition between stationary and periodic behavior. If 
the system is in a stable equilibrium before the bifurcation, the stability is lost at the bifurcation point. As a result 
the system abundance of species start to oscillate periodically. And this situation occurs only when the diffusion 
vanishes. Mathematically speaking, the Hopf bifurcation occurs when Im(A(fc)) ^ 0, Re(A(fc)) = at k — 0.. Then 
we can get the critical value of the transition, Hopf bifurcation parameter — K, equals 

K = B (we -g - I3e-rj) 2 

(we — 1) (r; 2 — (3 2 e 2 + rwe 2 (3 + e 2 (3 w-q — werj 2 ) 

At the Hopf bifurcation threshold, the temporal symmetry of the system is broken and gives rise to uniform 
oscillations in space and periodic oscillations in time with the frequency 

ujh = Im(A(/c)) = v'Ao, 
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where 

I K(r)-P e)(K(rwe-0 e+rj) 2 +7)6-0 8 e+A rwe rj B+r8 e w)+2 r5 r\ Be w 1 rj 

A _ __V J_ 

^° 0Ke 2 w(rKwe-f3Ke+Krj+S) ' 

and 

/ \ V 2 

S = [r 2 K 2 w 2 e 2 - 2 rK 2 ws 2 [3 + 2 rK 2 we 77 + f3 2 K 2 e 2 -2(3K 2 er] + K 2 rf + 4 rwe Kr, Bj . 
The corresponding wavelength is 

. 2tt 2tt 
\ H = 



The Turing instability is not dependent upon the geometry of the system but only upon the reaction rates and 
diffusion. And it cannot be expected when the diffusion term is absent and it can occur only when the activator (e.g., 
N) diffuses more slowly than the inhibitor (e.g., P). Linear analysis above shows that the necessary conditions for 
yielding Turing patterns are given by 

fn+9 P < 0, f n 9p - f P 9n > 0, d 2 g p + dif n > 0, (d 2 f n + dig p ) 2 > Ad\d 2 (j n g v - f P g n )- 

Mathematically speaking, as d\ <C d 2 , the Turing bifurcation occurs when 

Im(A(fc)) = 0, Re(A(fc)) = at k = k T ^ 0, 



and the wavenumber kr satisfies k\ = y ■ At the Turing bifurcation threshold, the spatial symmetry of the system 
is broken and the patterns are stationary in time and oscillatory in space with the corresponding wavelength 

A T = g. (8) 



And the critical value of Turing bifurcation parameter K takes the following form: 

FiA+F 2 
G1A+G2 

where 



Kt = (9) 



M 2 



A = (—r/di (ef3 — rj){—£ d 2 /3 — rjd 2 + ewrjdi) (— er] d\rw — erj d\ (3 + e r 2 d 2 w — e rd 2 f3 + d\ rf — rd 2 rf)) 

Fi = - [r) 2 {diew - d 2 ) 2 + sd 2 (3{ed 2 [3 + 2r)d 2 + Sew^dx)) 2 ((-2e 2 (3 wd 2 (2 wsifd 1 d 2 - Awe 2 rj d x d 2 {3 - rj 2 d 1 2 e 2 w 2 
+(3 2 d 2 2 e 2 - d 2 2 r, 2 )r) + 2 e d 2 (3 (e (3 - rf) (e 2 (d 2 (3 + w V d x ) 2 + 77 d 2 (rj d 2 - 2 e wr, d 1 + 2 e d 2 (3)))rB, 

F 2 = -(rj 2 (diew - d 2 ) 2 + ed 2 (3 (e d 2 (3 + 2r]d 2 + 6 e wi] di)) 2 (77 (3 e 2 wd 2 (e d 2 (3 + T]d 2 - e wq d 1 )(e 3 diw(3 d 2 (3 + wqdi) 2 
+d 2 3 (r] + e(3) 2 - e dir] 2 wd 2 (die w + d 2 ))r 2 - 77 (e (3 ~ i])(sd 2 (3 + rjd 2 - ewr) di)(e 4 i] 3 w 4 d 1 4: + e 3 rfw 3 d 2 {-Ari 
+9 e /3)di 3 + 3 e 2 r, w 2 d 2 2 (~5 77 e (3 + 2 rj 2 + 9 e 2 l3 2 )d x 2 + e wd 2 3 (ll e (3 - 4 77) (77 + e f3) 2 di + d 2 4 (r] + e f3) 3 )r 
+2i 1 f3d l ed 2 {ef3-vi) 2 {ed 2 f3 + r l d 2 -ewr l d 1 ){e 2 {d 2 l3 + wr 1 d 1 ) 2 +r l d 2 {r l d 2 - 2ewr]d 1 + 2ed 2 f3)))rB, 

G\ = G11G12, 

Gn = <ied 2 (3(-(£ 2 (3d 2 (d 2 (3 -Awr)di)- r) 2 (diew- d 2 ) 2 )wer + {e (3 - rj){e 2 {d 2 (3 + wr/di) 2 + r\ d 2 (rj d 2 - 2 e w-q dx 
+2ed 2 /3))), 

G12 = (e 3 diw(3d 2 f3 + wr / d 1 ) 2 + d 2 3 (77 + e f3) 2 - e diTj 2 wd 2 (die w + d 2 ))d 2 w/3 e 2 r 2 - (e (3 - rj)(e 3 d 2 3 (d 2 + 11 d x e w)[3 3 
+3 d 2 T) e 2 (3 die w + d 2 ) 2 f3 2 + 3 e r/ 2 d 2 (3 d x e w + d 2 )(die w - d 2 ) 2 (3 + rf^e w - d 2 ) 4 )r + 2 e (3 d^e (3 - t/) 2 
{e 2 {d 2 (3 + wr] d{) 2 + 77 d 2 (77 d 2 - 2 e wi] d 1 + 2 e d 2 (3)) , 

G 2 = G21G22, 

G 2 i = ed 2 f3 + rjd 2 - ewrjdx, 

G 22 = 3o + gir + g 2 r 2 + g 3 r 3 + g 4 r 4 , 

go = 8 e 2 /3 2 77 d^d^ie (3 - r/) 4 (e 2 (d 2 /3 + wr 1 d 1 ) 2 -2e d 1 r 1 2 wd 2 + d 2 2 r/ 2 + 2 (3 e 77 d 2 2 ) 2 , 

g! = Ae (3d 1 d 2 {e (3 - r]) 3 {e 2 {d 2 (3 + wrjdx) 2 -2ed 1 rj 2 wd 2 + d 2 r) 2 + 2 (3 e 77 d 2 2 )(-(d 1 e w - d 2 frf - 2ej3d 2 {Zd 1 sw + d 2 ) 
(d l£ w- d 2 ) 2 v 3 - 16 e 3 f3 2 d 1 wd 2 2 (d 1 ew + d 2 )rj 2 - 2e 3 d 2 3 (3 3 {-d 2 + 5d l£ w)ri + d 2 4 e 4 f3 4 ), 
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.92 = (e - v) 2 ((diew - d 2 ) V + 6 e d 2 (3 d x e w + d 2 ){d 1 e w - d 2 frf + 2 d 2 2 e 2 (151 e 2 w 2 d 1 2 + 106 e diwd 2 + 15 d 2 ) 
{d x ew- d 2 ) 4 r) 5 + 4e 3 d 2 3 3 (156s 3 d 1 3 w 3 + I33s 2 w 2 d 1 2 d 2 + ^wd x ed 2 2 + 5d 2 3 ){d 1 sw ~ d 2 ) 2 rf + d 2 4 s 4 4 {3d 1 ew 
+d 2 )(389e 3 d 1 3 w 3 + 117 e 2 w 2 d 1 2 d 2 + 183wdied 2 2 + 15d 2 3 )ry 3 + 2 e 5 /3 5 d 2 5 (357 e 3 d 1 3 w 3 + 297 s 2 w 2 d 1 2 d 2 
+47 wd l£ d 2 2 + 3 d 2 3 )r; 2 + e 6 /3 6 d 2 6 (153 s 2 w 2 d 1 2 + d 2 - 10 e d lW d 2 )r] - 12 e 8 P 7 wd 1 d 2 7 ), 

,g 3 = -2e 2 l3wd 2 (e (3 - 77) ((die w + d 2 )(c?ieu; - d 2 )V + e/3d 2 (13e 2 w 2 c?i 2 + Qediwd 2 + 5d 2 2 )(d 1 ew - d 2 ) 4 ?y 5 
+2/3 2 d 2 2 e 2 (38e 3 di 3 w 3 + e 2 w 2 di 2 d 2 + 20wdied 2 2 + 5d 2 3 )(d l£ w - d 2 )V + 2e 3 /3 3 rf 2 3 (5d 2 4 + 117e 4 w 4 di 4 
+Ued lW d 2 3 + Ud 2 2 d 1 2 w 2 e 2 - 52 d 2 d 1 3 w 3 e 3 )rf + d 2 A e 4 ft 4 {227 e 2 w 2 d 1 2 d 2 + 329 e 3 d 1 3 w 3 + 19wd Y ed 2 
+5d 2 3 )r] 2 + s 5 (3 5 d 2 5 (Ued 1 wd 2 + d 2 + 121 s 2 w 2 d 1 2 )r 1 - 6 e 7 d lW d 2 6 f3 6 ) , 

.94 = e 4 /3 2 w 2 rf 2 2 ((e 2 w 2 di 2 + 6ed lW d 2 + d 2 2 )(wd l£ - d 2 ) 4 r/ 5 + As f3d 2 {2e 3 d 1 3 w 3 + lie 2 w 2 d 2 d 2 + d 2 3 )(wd l£ - d 2 ) 2 r, 4 
+2fJ 2 d 2 2 e 2 (3d 2 4 + 52d 2 d! 3 w 3 e 3 - 6d 2 2 d 1 2 w 2 e 2 + Aed lW d 2 3 + lle^d^ 3 + Ae 3 (3 3 d 2 3 {-9e 2 w 2 d 1 2 d 2 + d 2 3 
+ 13 e 3 d 1 3 w 3 + 11 wdxe d 2 2 )r] 2 + d 2 V/3 4 (113 e 2 w 2 di 2 + d 2 2 + 22 e d lW d 2 )r] - 4 d 2 5 s 6 wd 1 (3 5 ), 

In the following, linear stability analysis yields the bifurcation diagram with r = 0.5, e = 1, (3 = 0.6, B = 0.4, 
T] = 0.25, w = 0A 7 d 2 = 1 shown in Fig. EUA). 

The Hopf bifurcation line and the Turing bifurcation curve separate the parametric space into four distinct domains. 
In domain I, located below all two bifurcation lines, the steady state is the only stable solution of the system. Domain 
II is the region of pure Turing instability, while domain III is the region of pure Hopf instability. In domain IV, which 
is located above all two bifurcation lines, both Hopf and Turing instability occur. 

To see the relation between the real and the imaginary parts of the eigenvalue X(k), we plot in Fig. HJB)-(F) the 
real and the imaginary parts of the eigenvalue at different K with r = 0.5, e = 1, /3 = 0.6, B = 0.4, 77 = 0.25, w = 0.4 
di = 0.01 and d 2 = 1 for the system [5] 




FIG. 1: (A) K — d\ Bifurcation diagram for the model[2]with r = 0.5, e = 1, /3 = 0.6, B = 0.4, 77 = 0.25, w = 0.4 and d 2 = 1. 
Hopf and Turing bifurcation lines separate the parameter space into four domains. And the Hopf- Turing bifurcation point is 
(di,K) — (0.02742,2.63158). In (B) — (F), Re(A(fe)) and Im(A(fe)) are shown by solid curves and dotted curves, respectively. 
The other parameters are: di = 0.01, the bifurcation parameter K: (B) 2.131712170, the critical value of Kt); (C) 2.2; (D) 
2.6; (E) 2.631578947 the critical value of K H ); (F) 3.0. 

From the definition of Hopf and Turing bifurcation, we know that the relation between the real, the imaginary parts 
of the eigenvalue A(fc) determine the bifurcation type. The relation between Re(A(/c)), Im(A(fc)) and k are shown in 
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figure QJB)-(F). Figure HJB) illustrate the case of parameter locate in domain I in figure QJA), K = 2.131712170, 
the critical value of Turing bifurcation, in this case, Re(A(fc)) > at k = while Im(A(fc)) ^ 0. In figure []JC)(D), 
K = 2.2 and K = 2.6, the parameter locate in domain II, the pure Turing instability occurs, one can see that at 
k = 0, Re(A(fc)) = 0, Im(A(fc)) ^ 0. Figure HfE), K = 2.631578947, the critical value of Hopf bifurcation, in this case, 
Re(A(fc)) = at k — while Im(A(fc)) ^ 0. When K = 3.0, parameter locate in domain IV, figure QjF) indicate that 
at k = 0, Re(A(fc)) > 0, Im(A(fc)) ^ 0. 



3. SPATIOTEMPORAL PATTERN FORMATION 



In this section, we perform extensive numerical simulations of the spatially extended model Q in two-dimensional 
spaces, and the qualitative results are shown here. All our numerical simulations employ the periodic Neumann 
(zero-flux) boundary conditions with a system size of 200 x 200 space units and r = 0.5, e — 1, (3 = 0.6, B — 0.4, 
X] = 0.25, w — 0.4 d\ = 0.01 and d,2 = 1. The EqJ5] are solved numerically in two-dimensional space using a finite 
difference approximation for the spatial derivatives and an explicit Euler method for the time integration with a time 
stepsize of At = 0.01 and space stepsize (lattice constant) Ah = 0.25 (see, for details, [lj|). The scale of the space 
and time are average to the Euler method. The initial density distribution corresponds to random perturbations 
around the stationary state (N* , P*) in model ([2]) with a variance significantly lower than the amplitude of the final 
patterns, which seems to be more general from the biological point of view. When the system reached a stable state 
(stationary or oscillatory), we took a snapshot with yellow levels linearly proportional to the free species density and 
red corresponding to high while blue corresponding to low. 

In the numerical simulations, different types of dynamics are observed and we have found that the distributions of 
predator and prey are always of the same type. Consequently, we can restrict our analysis of pattern formation to 
one distribution. In this section, we show the distribution of prey, for instance. 




. Vj : • : • * 




• • • »• • • • • 



FIG. 2: Dynamics of the time evolution of the prey of model [2] with di = 0.01, Kt < K = 2.2 < Kh- (A) iteration, (B) 
10000 iterations, (C) 200000 iterations, (D) 300000 iterations. 



Fig. [5] shows the evolution of the spatial pattern of prey at 0, 10000, 200000 and 300000 iterations, with random 
small perturbation of the stationary solution (N* , P*) — (0.5094, 0.7829) of the spatially homogeneous systems when 
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FIG. 3: Dynamics of the time evolution of the prey of model [2] with d\ = 0.01, Kt < K = 2.6 < Kh- (A) iteration, (B) 
10000 iterations, (C) 30000 iterations, (D) 100000 iterations. 




FIG. 4: Dynamics of the time evolution of the prey of model [5] with d\ = 0.01, Kt < Kh < K — 3.0. (A) iteration, (B) 
20000 iterations, (C) 50000 iterations, (D) 100000 iterations. 
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FIG. 5: Dynamics of the time evolution of the prey of model [2] with d\ = 0.04, Kt < Kh < K = 2.65. (A) iteration, (B) 
30000 iterations, (C) 70000 iterations, (D) 100000 iterations. 




FIG. 6: Dynamical behaviors of model [2] (A) Phase portrait; (B) Time-series plot; (C) Space-time plots corresponding to 
Fig. [5jC). The parameters are the same as those in Fig [5] 

K = 2.2, located in domain II, slightly more than the Turing bifurcation threshold Kt = 2.1317 and less than the 
Hopf bifurcation threshold Kh = 2.6316. In this case, one can see that for models the random initial distribution 
(c.f., Fig.fS^A)) leads to the formation of a regular macroscopic spotted pattern which prevails over the whole domain 
at last, and the dynamics of the system does not undergo any further changes (c.f., Fig.^D)). 

Fig. [3] shows the evolution of the spatial pattern of prey at 0, 10000, 30000 and 100000 iterations when K = 2.6 
which is more than Kt and slightly less than Kh- Although the dynamics of the system starts from the stationary 
solution (N*,P*) = (0.5252,0.8382), there is an essential difference in the case above. From the snapshots, one can 
see that the steady state of spotted pattern and the stripe- like pattern coexist (c.f., Fig. [3JD)). 

In Fig. 31 Kh < K — 3.0, i.e., parameters in domain IV, both Hopf and Turing instability occur. In this case, 
(N*,P*) = (0.5380,0.8831). One can see that the evolution of the spatial pattern of prey at 0, 20000, 50000 and 
100000 iterations. After the spatial chaos patterns (c.f., Fig. IDJB)), a regular stationary stripe-like spatial state 
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emerges (c.f., Fig.QjD)). 

When K = 2.65, d\ — 0.04, i.e., parameters in domain III (c.f., Fig. QJA)), pure Hopf instability occurs. As an 
example, the formation of a regular macroscopic two-dimensional spatial pattern, the spiral pattern, is shown in 
Fig. [5] with a system size of 400 x 400 space units. One can see that for model [H the random initial distribution 
around the steady state (N* , P*) = (0.5270, 0.8443), a unstable focus of model [21 leads to the formation of the spiral 
pattern in the domain (c.f., Fig.^D)). In other words, in this situation, spatially uniform steady-state predator-prey 
coexistence is no longer. Small random fluctuations will be strongly amplified by diffusion, leading to nonuniform 
population distributions. From the analysis in section II, we find with these parameters in domain III, the pattern 
formation, i.e., the spiral pattern, arises from pure Hopf instability. In order to make it clearer, in Fig. [SJ we show 
phase portrait (Fig. [S^A)), time series plot (Fig. OB)) and space-time plots (Fig. IUJC)), a one-dimensional example 
corresponding to Fig.^C). And the method of space-time plots is that let y be a constant (here, y = 200, the center 
line of each snapshots), from each pattern snapshots, choose the line y = 200, and pile these lines in-time-order. The 
space-time plots show the evolution process of the prey N throughout time t and space x. Fig. [6jA) exhibits the 
"local" phase plane of the system obtained in a fixed point (0.5270, 0.8443) inside the region invaded by the irregular 
spatiotemporal oscillations, and the trajectory fills nearly the whole domain inside the limit cycle. Fig.[6jB) illustrates 
the evolution process of prey density with time, periodic oscillations in time finally. From Fig. (6]JC), one can clearly 
see that a spiral wave emerges, and the system gives rise to uniform oscillations in space and periodic oscillations in 
time. 

Comparing Fig. [2j Fig. [3] with Fig. [4] Fig. [5j we can see that the bifurcation parameter K determines the type of 
the pattern formation even with the same parameters, e.g., r, e, /3, B, rj, w, d,2- In domain II of Fig. HJA), the closer 
K is to Kt, the more distinct the spotted spatial pattern becomes (c.f., Fig. [2^D)). When K is much closer to Kh, 
the spotted and stripe-like patterns coexist (c.f., Fig.^D)). When K is bigger than Kh and far away from the Turing 
bifurcation value Kt, the distinct stripe-like pattern emerges. When K is bigger than Kh and smaller than Kt, the 
spiral wave pattern occurs (c.f., Fig. [5jC,D)). So we may draw a conclusion that for model [2] pure Turing instability 
gives birth to the spotted pattern, pure Hopf bifurcation gives birth to the spiral wave pattern, and both of them give 
birth to the stripe-like pattern. 

4. CONCLUSIONS AND REMARKS 

In this paper, we have presented a theoretical analysis of evolutionary processes that involves organisms distribution 
and their interaction of spatially distributed population with local diffusion. And the numerical simulations were 
consistent with the predictions drawn from the bifurcation analysis, i.e., Hopf bifurcation and Turing instability. In 
the domain II of Fig. [TJ A) , the stationary state of periodic spotted pattern exists when the parameters are near the 
Turing bifurcation line, while near the Hopf bifurcation line, both the spotted pattern and the stripe-like pattern 
coexist. When the parameters are located in domain III, pure Hopf instability occurs, and the spiral wave pattern 
emerges. When the parameters are located in domain IV, both Hopf and Turing instability occur, and the stationary 
state of stripe-like pattern exists. 

Turing instability is relevant not only in reaction-diffusion systems, but also in describing other dissipative structures, 
which can be understood in terms of diffusion-driven instability. In addition to the biological relevance of Turing 
systems, their ability to generate structure is of great interest from the point of view of physics. There are various 
physical systems that show similar phenomena, although the underlying mechanisms can be very different. Thus, most 
of the research in the field relies on experiments and numerical simulations justified by an analytical examination. In 
addition, in simulations one may study pattern formation under constraints that are beyond the reach of experiments 
and the numerical data is also easy to analyze. 

In Ref. 0, it's indicated that the basic idea of diffusion-driven instability in a reaction-diffusion model can be 
understood in terms of an activator-inhibitor system. And a random increase of activator species (prey, N) has a 
positive effect on the creation rate of both activator and inhibitor (predator, P) species. In other words, random 
fluctuations may cause a nonuniform prey density. This elevated prey density has a positive effect both on prey and 
predator population growth rates. Following D. Alonso et al Q, we give the discussion to model([2]). From Eqs.([2]), 
we can obtain the following equations: 

1 dN _ /-t N\ PP 1 8P _ ef3N nm 

N dt 'V- 1 - K> B+N+wPi P dt B+N+wP '/• V ±u ) 

Similar to Ref. the first equation in Eqs. (|10[) is a one- humped function of prey density, the numerical result can 
be found in (20j , and prey growth rate can be increased by a higher local prey density at least in a range of parameter 
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values. On the other hand, the second equation in Eqs. (|10[) . i.e., predator numerical response, is an ever-increasing 
function of N, and high prey density always has a positive influence on predator growth. More importantly, inhibitor 
species (predator, P) must diffuse faster than activator species (prey, N), for an increment in inhibitor species may 
have a negative effect on formation rate of both species. Thus, as random fluctuations increase local prey density 
over its equilibrium value, prey population undergoes an accelerated growth. Simultaneously, predator population 
also increases, but as predators diffuse faster than prey, they disperse away from the center of prey outbreaks. If 
relative diffusion (d^/di) is large enough, prey growth rate will reach negative values and prey population will be 
driven by predators to a very low level in those regions. The final result is the formation of patches of high prey 
density surrounded by areas of low prey density. Predators follow the same pattern. 

In Ref. [ll|, Neuhauser and Pacala formulated the Lotka-Volterra model as a spatial one. They found the striking 
result that the coexistence of patterns is actually harder to get in the spatial model than in the non-spatial one. One 
reason can be traced to how local interaction between individual members of the species are represented in the model. 
In this study, our results show that the predator-prey model with Beddington-DeAngelis-type functional response 
and reaction-diffusion (e.g., Eqs.EJ also represents rich spatial dynamics, such as spotted pattern, stripe-like pattern, 
coexistence of both stripe-like and spotted pattern, spiral pattern, etc. It will be useful for studying the dynamic 
complexity of ecosystems or physical systems. In Ref. 9], it is indicated that reaction-diffusion models provide a way 
to translate local assumptions about the movement, mortality, and reproduction of individuals into global conclusions 
about the persistence or extinction of populations and the coexistence of interacting species. They can be derived 
mechanistically via rescaling from models of individual movement which are based on random walks, i.e., small random 
perturbation of the stationary solution (N*,P*) of the spatially homogeneous model j2]). Reaction-diffusion system, 
i.e., model ([2]), is spatially explicit and typically incorporate quantities such as dispersal rates, local growth rates, and 
carrying capacities as parameters which may vary with location or time. 

More interesting, when the parameters are located in domain III, pure Hopf instability occurs, and the spiral wave 
pattern emerges (c.f., Fig. [5l Fig. [6^0)). To our knowledge, we haven't got any report about one system that has 
spotted, stripe-like and spiral pattern formation meantime. It's well known that spirals and curves are the most 
fascinating clusters to emerge from the predator-prey model. A spiral will form from a wave front when the rabbit 
line (which is leading the front) overlaps the pursuing line of predator. The prey on the extreme end of the line stop 
moving as there are no predator in their immediate vicinity. However the prey and the predator in the center of the 
line continue moving forward. This forms a small trail of prey at one (or both) ends of the front. These prey start 
breeding and the trailing line of prey thickens and attracts the attention of predator at the end of the fox line that 
turn towards this new source of prey. Thus a spiral forms with predator on the inside and prey on the outside. If the 
original overlap of prey occurs at both ends of the line a double spiral will form. Spirals can also form as a prey blob 



collapses after predator eat into it [23. |24| . Thus, reaction-diffusion system provides a good framework for studying 
questions about the ways that habitat geometry and the size or variation in vital parameters influence population 
dynamics. 

On the other hand, one can see that although there is a small difference between the denominators of the functional 
responses of Michaelis-Menten-type and those of Beddington-DeAngelis-type, there is an enormous gap between them 
in the process of computations. The Turing bifurcation expression of Michaelis-Menten-type predator-prey system is 
simple [641 ], while from @, one can see that the Turing bifurcation analysis requires huge-sized computations, so we 
have to obtain more help via computers. 

In fact, computer aided analysis is useful for nonlinear analysis. And computers have played an important role 
throughout the history of ecology. Today, numerical simulations also play an important role in spatial ecology. There 
are some international mathematical softwares, such as Matlab, Maple, Mathematica, etc. We have finished all our 
symbolic computations in Maple and obtained our pattern snapshots (i.e., numerical simulations) in Matlab for Maple 
is more superior in symbolic computations while Matlab is more superior in numerical computations. 
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In this paper, we investigate the emergence of a predator-prey model with Beddington-DeAngelis- 
type functional response and reaction-diffusion. We derive the conditions for Hopf and Turing 
bifurcation on the spatial domain. Based on the stability and bifurcation analysis, we give the 
, spatial pattern formation via numerical simulation, i.e., the evolution process of the model near the 

coexistence equilibrium point. We find that for the model we consider, pure Turing instability gives 
1^-^ ' birth to the spotted pattern, pure Hopf instability gives birth to the spiral wave pattern, and both 

Hopf and Turing instability give birth to stripe-like pattern. Our results show that reaction-diffusion 
i i model is an appropriate tool for investigating fundamental mechanism of complex spatiotemporal 

PQ ' dynamics. It will be useful for studying the dynamic complexity of ecosystems. 
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Si . 1. INTRODUCTION 
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Mathematical models have played an important role throughout the history of ecology. Models in ecology serve 
a variety of purposes, which range from illustrating an idea to parameterizing a complex real-world situation. They 
are used to make general predictions, to guide management practices, and to provide a basis for the development of 
statistical tools and testable hypotheses (ToL 4fl 51 1 . 
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A fundamental goal of theoretical ecology is to understand how the interactions of individual organisms with each 
other and with the environment determine the distribution of populations and the structure of communities 0] . As 
we know, our ecological environment is a huge and highly complex system. This complexity arises in part from 



the diversity of biological species, and also from the complexity of every individual organism 2^, 3^]. The dynamic 
behavior of predator-prey model has long been and will continue to be one of the dominant themes in both ecology 
and mathematical ecology due to its universal existence and importance 0, • 

Before the 1970s, ecological population models typically used ordinary differential equations, seeking equilibria and 
analyzing stability. The early models provided important insights, such as when species can stably coexist and when 
predator and prey density oscillate over time [loj ]. 
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We live in a spatial world, and spatial patterns are ubiquitous in nature, which modify the temporal dynamics and 
stability properties of population density at a range of spatial scales, whose effects must be incorporated in temporal 
ecological models that do not represent space explicitly. And the spatial component of ecological interactions has 
been identified as an important factor in how ecological communities are shaped (Io| . Empirical evidence suggests 
that the spatial scale and structure of the environment can influence population interactions and the composition of 
communities [§]. In recent decades, the role of spatial effects in maintaining biodiversity has received a great deal of 

attention in the literature on conservation d m m, h 0, m m m . 

The past investigations have revealed that spatial inhomogcneities like the inhomogeneous distribution of nutrients 
as well as interactions on spatial scales like migration can have an important impact on the dynamics of ecological 
populations [Hoi . [Hlj . In particular it has been shown that spatial inhomogeneities promote the persistence of ecological 
populations, play an important role in speciation and stabilize population levels Spatial ecology today is still 
dominated by theoretical investigations, and empirical studies that explore the role of space are becoming more 
common due to technological advances that allow the recording of exact spatial locations [l0| • 

In 1952, Alan Turing, one of the key scientists of 20th century, mathematically showed that a system of coupled 
reaction-diffusion equations could give rise to spatial concentration patterns of a fixed characteristic length from an 
arbitrary initial configuration due to diffusion-driven instability [611 ]. The work by Turing belongs to the field of 
pattern formation, a subficld of mathematical biology. Pattern formation in nonlinear complex systems is one of 
the central problems of the natural, social, and technological sciences. The occurrence of multiple steady states and 
transitions from one to another after critical fluctuations, the phenomena of excitability, oscillations, waves and the 
emergence of macroscopic order from microscopic interactions in various nonlinear noncquilibrium systems in nature 
and society have been the subject of many theoretical and experimental studies [50j]. It has been qualitatively shown 
that Turing models can indeed imitate biological patterns [37| . 

The study of biological pattern formation has gained popularity since the 1970s, and Segel and Jackson [56j were 
the first to apply Turing's ideas to a problem in population dynamics: the dissipative instability in the predator-prey 
interaction of phytoplankton and herbivorous copepods with higher herbivore motility. At the same time, Gierer and 
Meinhardt [2 11 ] gave a biologically justified formulation of a Turing model and studied its properties by employing 
numerical simulations. Levin and Segel [HI ] suggested this scenario of spatial pattern formation was a possible origin 
of planktonic patchiness. 

Spatial patterns and aggregated population distributions are common in nature and in a variety of spatiotemporal 
models with local ecological interactions 54] . And the understanding of patterns and mechanism of species spatial 
dispersal is an issue of current interest in conservation biology and ecology [3, 39, 50]. It arises from many ecological 
applications, and in p articular, p lays a major role in connection to biological invasion and epidemic spread and 
so on 0, @, [n], 0, 33, 43, 47, 43, Ho, 51 1. The field of research on pattern formation modeled by reaction-diffusion 
systems, which provides a general theoretical framework for describi ng p attern formation in systems from many diverse 
disciplines includingbiology @, 0, IMIMJisl M, M, M, M, H], chemistry [13, M, M, M, M, E3, El, M, , 



physics [12, 131|, |38|, |60| , epidemiology [43 
during the last decade. 

In general, a classical predator-prey model can be written as the form 0, 0] 



and so on, seems to be a new increasingly interesting area, particularly 



N = Nf(N) - Pg(N,P), P = h\g{N,P),P]P, 

where N and P are prey and predator densities, respectively, f(N) the prey growth rate, g(N, P) the functional 
response, e.g., the prey consumption rate by an average single predator, and h[g(N, P), P] the per capita growth rate 
of predators (also known as the "predator numerical response"), which obviously increases with the prey consumption 
rate. The most widely accepted assumption for the numerical response is the linear one Q : 

h\g(N,P),P] = eg(N,P)-ri 



where r) is a per capita predator death rate and e the conversion efficiency of food into offspring. 

In population dynamics, a functional response g(N, P) of the predator to the prey density refers to the change in 
the density of prey attached per unit time per predator as the prey density changes [jl [Hj ■ There have been several 
famous functional response types: Holling types I — III 0, [S3]; Hassell-Varley type Eqk Beddington-DeAngelis type 
by Bcddington [6j and DeAngelis et al 15] independently; the Crowley-Martin type [l3j; and the recent well-known 
ratio-dependence type by Arditi and Ginzburg Q later studied by Kuang and Beretta [35j . Of them, the Holling 
type I-III is labeled "prey-dependent" and other types that consider the interference among predators are labeled 
"predator-dependent" 0] . 
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In our previous work [64| , we studied the spatiotemporal complexity of a ratio-dependent predator-prey model with 
Michaelis-Menten-type functional response. Compare Michaelis-Menten-type functional response 

(a, h are positive constants, a capture rate and h handling time) with Bcddington-DeAngelis-type functional response 

g (N, P) = ^ (1) 

JK ' ' B + N + wP w 

(/?, B are positive constants, (3 a maximum consumption rate, B a saturation constant, w a predator interference 
parameter, a constant, w < is the case where predators benefit from cofeeding). Some scholars 0, 46|, 5t| indicat 



that where there is a small difference between the denominators, there is a world between them from the biological 
point. 

If predators do not waste time interacting with one another or if their attacks are always successful and instanta- 
neous, i.e., w — in ( [1} , then a Holling type II functional response is obtained: 

m- 0N 



B + N 



If B = in (TT|), the Michaelis-Menten-type functional response is obtained. 

Furthermore, in (59j . by comparing the statistical evidence from nineteen categories of predator- prey systems with 
three predator-dependent functional responses, Skalski and Gilliam pointed out that the predator-dependent functional 
responses could provide better descriptions of predator feeding over a range of predator-prey abundance, and in some 
cases, the Beddington-DeAngelis-type functional response performed even better [Zj|. 

Some progress has been seen in the study of predator-prey model with Beddington-DeAngelis-type functional 
response 0, |, Q [H, H [H HI, 0, [H |H IS, El . However, the research on considering reaction-diffusion to such 



a model, to our knowledge, seems rare. 

In this paper, we report a study of Turing pattern formation in a two-species reaction-diffusion predator-prey model 
with Beddington-DeAngelis-type functional response. In the next section we give a brief stability and bifurcation 
analysis of the model. Then, we present and discuss the results of numerical simulations, which is followed by the last 
section, i.e., conclusions and remarks. 



2. STABILITY AND BIFURCATION ANALYSIS 



In this paper, we mainly focus on the following predator- prey model with Beddington-DeAngelis-type functional 
response and reaction-diffusion: 

^ = r (1 - f ) N - ^§L_P + d^N, % = -r,P + d 2 ^P (2) 



where t denotes time and N, P stand for prey and predator density, respectively. All parameters are positive constants, 
r standing for maximum per capita growth rate of the prey, j3 capture rate, r\ predator death rate, w a predator 
interference parameter and K carrying capacity, which is the nonzero equilibrium population size. The diffusion 
coefficients are denoted by d\ and c?2, respectively. V 2 = + -gjp is the usual Laplacian operator in two-dimensional 
space. 

The first step in analyzing the model is to determine the equilibria (stationary states) of the non-spatial model 
obtained by setting space derivatives equal to zero, i.e., 

1 V 1 A'/ iv B+N+wP r u > B+N+wP r 'l r u - y°> 

In fact, physically, an equilibrium represents a situation without "life". It may mean no motion of a pendulum, no 
reaction in a reactor, no nerve activity, no flutter of an airfoil, no laser operation, or no circadian rhythms of biological 
clocks [57j | . And at each equilibrium point, the movement of the population dynamics vanishes. 

Eqs.® has at most three equilibria, which correspond to spatially homogeneous equilibria of the full model, Eq.([2]), 
in the positive quadrant: 

(i) (0,0) (total extinct) is a saddle point. 
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(ii) (K, 0) (extinct of the predator, or prey-only) is a stable node if e j3 < rj or if e f3 > r\ and K < — _^p +v ', a saddle 
if e (3 < T) and K > — zf^p" ! a saddle- node if e j3 < rj and K — — ■ 

(hi) a nontrivial stationary state (N*,P*) (coexistence of prey and predator), where 

N* = — \K(rwe -e(3+ rj) + y/K 2 (rwe - e (3 + rf) 2 + 4 rKwerjB^ , 

p* _ (0e-y) jy* _ B 
wt] w ' 

To perform a linear stability analysis, we linearize the dynamic system [2] around the equilibrium point (N*, P*) for 
small space- and time-dependent fluctuations and expand them in Fourier space 

N(x, t) - N*e xt e lts , P(x, t) - P*e xt e^ s , 
and obtain the characteristic equation 

\A~k 2 D-\I\ =0, (4) 
where D = diag(di, di) and the Jacobian matrix A is given by 

9„f d p f \ ( f n fp \ 

(N*,P*) V 9n 9p J 

Now Eq.(|3]) can be solved, yielding the so called characteristic polynomial of the original problem 

A 2 -tr fe A + Afc = 0, (5) 

where 

tr fc = f n +g p - k 2 (dx + d 2 ) = tr - k 2 (d 1 + d 2 ), 

Afc = f„g p - f p g„ - k 2 (f n d 2 + g p dx) + k 4 did 2 = A - k 2 (f n d 2 + g p di) + k 4 dxd 2 . 
The roots of Eq. ([5|) yield the dispersion relation 



A = 

\d n g d p g 



A(fc) = ^(tr fe ± x /tr2-4A fc ). (6) 



We know that one type of instability (or bifurcation) will break one type of symmetry of a system, i.e., in the 
bifurcation point, two equilibrium states intersect and exchange their stability. Biologically speaking, this bifurcation 
corresponds to a smooth transition between equilibrium states. The Hopf bifurcation is space-independent and it 
breaks the temporal symmetry of a system and gives rise to oscillations that are uniform in space and periodic in 
time. The Turing bifurcation breaks spatial symmetry, leading to the formation of patterns that are stationary in 
time and oscillatory in space. 

The Hopf instability or bifurcation is an important instability in reaction-diffusion systems for which the conditions 
result in a stable limit cycle (oscillations). In terms of the linearized problem (EqHJ) the onset of Hopf instability 
corresponds to the case, when a pair of imaginary eigenvalues cross the real axis from the negative to the positive side. 
The Hopf bifurcation of a equilibrium state is reflected by a transition between stationary and periodic behavior. If 
the system is in a stable equilibrium before the bifurcation, the stability is lost at the bifurcation point. As a result 
the system abundance of species start to oscillate periodically. And this situation occurs only when the diffusion 
vanishes. Mathematically speaking, the Hopf bifurcation occurs when Im(A(fc)) ^ 0, Re(A(fc)) = at k — 0.. Then 
we can get the critical value of the transition, Hopf bifurcation parameter — K, equals 

R = Bjweri- I3e-rj) 2 

(we — 1) (r; 2 — (3 2 e 2 + rwe 2 (3 + e 2 (3 wrj — we rj 2 ) 

At the Hopf bifurcation threshold, the temporal symmetry of the system is broken and gives rise to uniform 
oscillations in space and periodic oscillations in time with the frequency 

ujh = Im(A(/c)) = v'Ao, 
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where 

^K(rj— [3 e)(K(rwe — f3 e+r/) 2 +?7 5—f3 5 £+4 rwe r/ B+r5 e w)+2 rS r/ Be w^jrj 
^° = $ Ke 2 w(rKwe-f3 Ke+K-q+5) ' 

and 



S = (r 2 K 2 w 2 e 2 - 2 rK 2 we 2 /3 + 2 rK 2 we rj + f3 2 K 2 e 2 -2f3K 2 er] + K 2 rf + 4 rwe Kr] s) ^ . 
The corresponding wavelength is 



2tt 2tt 

As 



The Turing instability is not dependent upon the geometry of the system but only upon the reaction rates and 
diffusion. And it cannot be expected when the diffusion term is absent and it can occur only when the activator (e.g., 
N) diffuses more slowly than the inhibitor (e.g., P). Linear analysis above shows that the necessary conditions for 
yielding Turing patterns are given by 

/™+.g P <0, f n g p - f p g n > 0, d 2 g p + dif n >0, (d 2 f n + dig p ) 2 > 4did 2 (f n g p - f P 9n)- 

Mathematically speaking, as d\ <C d 2 , the Turing bifurcation occurs when 

Im(A(fc)) = 0, Rc(A(fc)) = at k = k T ^ 0, 



and the wavenumber kr satisfies k\ = y -jf^ . At the Turing bifurcation threshold, the spatial symmetry of the system 
is broken and the patterns are stationary in time and oscillatory in space with the corresponding wavelength 

A T = f T (8, 



And the critical value of Turing bifurcation parameter K takes the following form: 

F X A+F 2 
G1A+G2 

where 



k t - q a+g? 1 ( y i 



A = (—rjdi (e(3 — rj)(—ed 2 (3 — r\ d 2 + ewrjdi) 2 (—erjdi rw — er\ d\ (3 + er 2 d 2 w — e rd 2 (3 + d\ rj 2 — rd 2 77)) 1 2 , 

Fi = - {r] 2 {disw - d 2 ) 2 + ed 2 (3{ed 2 f3 + 2r)d 2 + &ewqd 1 )) 2 ((-2e 2 (3 wd 2 (2 we i] 2 did 2 - Awe 2 rj did 2 (3 - rj 2 d 1 2 e 2 w 2 
+(3 2 d 2 2 e 2 - d 2 2 r, 2 )r) + 2 e d 2 (3 (e (3 - rf) (e 2 (d 2 /3 + w V d x ) 2 + r, d 2 (r, d 2 - 2 e wr, d 1 + 2 e d 2 (3)))rB, 

F 2 = -(rj 2 (diew - d 2 ) 2 + ed 2 [3 (e d 2 (3 + 2r)d 2 + 6 e wi] di)) 2 (rj (3 e 2 wd 2 (e d 2 (3 + r/d 2 - e wi] di)(e 3 diw(3 d 2 (3 + wqdif 
+d 2 a (rj + e (3) 2 - e dir] 2 wd 2 (d 1 ew + d 2 ))r 2 - rj (e j3 - rj)(ed 2 (3 + rjd 2 - ewr] di){e 4 rfw i di A + e 3 r] 2 w 3 d 2 (-4:T] 
+9 e f3)d 1 3 + 3 e 2 r) w 2 d 2 2 (~5 r] e (3 + 2 rr 2 + 9 e 2 (3 2 )di 2 + e wd 2 3 (l 1 e (3 - 4 r?)(r ? + e (3) 2 d 1 + d 2 4 ( v + e f3) 3 )r 
+2i 1 f3d l ed 2 (e [3 - rj) 2 (e d 2 /3 + rj d 2 - ewqdi){e 2 {d 2 f3 + wr]di) 2 +r,d 2 (r/d 2 - 2 e wr) d x + 2 e d 2 (3)))rB, 

G\ = G11G12, 

Gn = 4ed 2 (3 (-(e 2 f3d 2 (d 2 f3 - Awr/di) - r) 2 (diew - d 2 ) 2 )wer + (e (3 - r))(e 2 (d 2 /3 + wrjdi) 2 + r\ d 2 (rj d 2 - 2 e wrj di 
+2ed 2 f3))), 

G12 = (e 3 diw(3d 2 /3 + wrjdi) 2 + d 2 3 (rj + e (3) 2 - e dirj 2 wd 2 (die w + d 2 ))d 2 wf3 e 2 r 2 - (e (3 - r])(e 3 d 2 3 (d 2 + 11 w)/3 3 
+3 d 2 2 r] e 2 (3 die w + d 2 ) 2 f3 2 + 3 e r) 2 d 2 (3 die w + d 2 ){die w - d 2 ) 2 (3 + rj 3 (die w - d 2 ) 4 )r + 2 e [3 did 2 (e (3 - r]) 2 
(e 2 (d 2 /3 + wr/ di) 2 + r) d 2 (r] d 2 - 2 e wr] di + 2 e d 2 (3)) , 
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^2 — ( - T 21 ( - T 22; 

G21 = e d 2 [3 + 77 d 2 — e wq di , 

G22 = 9o + 9if + 92 r 2 + g-3r 3 + 9^, 

.90 - 8 e 2 f3 2 i 1 d 1 2 d 2 2 (e (3 - r 1 ) i {e 2 {d 2 f3 + wi 1 d 1 ) 2 -2e d 1 rj 2 wd 2 + d 2 2 rf + 2 [3 e i] d 2 2 ) 2 , 

gi = 4e/?did 2 (£/3- ri) 3 (e 2 (d 2 /3 + wrjdx) 2 -2ed 1 ri 2 wd 2 + d 2 V + 2 f3 £ rj d 2 2 ){-{d 1 e w - d 2 ) i ri i - 2 e [3 d 2 (3 die w + d 2 ) 
{die w - d 2 ) 2 rf - 16 e 3 f3 2 d 1 wd 2 2 {d 1 £ w + d 2 )r] 2 - 2 e 3 d 2 3 (3 3 {-d 2 + 5 d x e w)rj + d 2 i e i (3 i ), 

g 2 = {ef3-r 1 ) 2 {{d 1 ew- d 2 )W + 6 e (3 d 2 {3 d x e w + d 2 ){d 1 e w - d 2 ) e r] 6 + f3 2 d 2 2 e 2 {151 e 2 w 2 d 2 + 106 e d x wd 2 + 15 d 2 2 ) 
{d l£ w~ d 2 ) 4 r] 5 + 4: e 3 d 2 3 p 3 {156 s 3 d 1 3 w 3 + 133£ 2 w 2 d 1 2 d 2 + 58 wd 1 sd 2 2 + 5d 2 3 ){d 1 sw - d 2 ) 2 rj i + d 2 4 s 4 [3 A (3 d x e w 
+d 2 ){389e 3 d 1 3 w 3 + 117 e 2 w 2 d 1 2 d 2 + 183wdied 2 2 + 15d 2 3 W + 2 e 5 f3 5 d 2 5 {357 e 3 d 1 3 w 3 + 297 e 2 w 2 d 1 2 d 2 
+47 wd l£ d 2 2 + 3 d 2 3 )ri 2 + e 6 /3 6 d 2 6 (153 s 2 w 2 d 1 2 + d 2 2 - 10 e d lW d 2 )r] - 12 e 8 f3 7 wd 1 d 2 7 ), 

g 3 = -2e 2 Pwd 2 {sP-r]){{d 1 ew + d 2 )(d 1 ew- d 2 frf + e(3d 2 {lZe 2 w 2 d 2 + §ediwd 2 + 5d 2 2 )(d 1 ew - d 2 ) 4 ry 5 
+2(3 2 d 2 2 e 2 {i8e 3 d 1 3 w 3 + e 2 w 2 d 1 2 d 2 + 2Qwd 1 ed 2 2 + 5d 2 3 ){d 1 ew~ d 2 ) 2 r 1 A + 2 e 3 (3 3 d 2 3 {5 d 2 A + 117e 4 w 4 di 4 
+44ediwd 2 3 + Ud 2 2 d 1 2 w 2 e 2 - 52 d 2 d x 3 w 3 e 3 )rf + d 2 4 e 4 ^ {227 s 2 w 2 d 1 2 d 2 + 329 e 3 d 1 3 w 3 + 79 wd x ed 2 2 
+5d 2 3 )r] 2 + e 5 f3 5 d 2 5 {Ued!wd 2 + d 2 2 + 121 e 2 w 2 d 1 2 )r ] - 6e 7 d lW d 2 6 l3 6 ), 

,g 4 = e 4 f3 2 w 2 d 2 2 {{e 2 w 2 d 1 2 + 6ed lW d 2 + d 2 2 ){wd l£ ~ d 2 )V + Ae P d 2 {2 e 3 d 1 3 w 3 + 13e 2 w 2 d 1 2 d 2 + d 2 3 ){wd l£ - d 2 )V 
+2/3 2 d 2 V(3d 2 4 + 52d 2 di 3 w 3 e 3 -6d 2 2 d 1 2 w 2 e 2 + 4ed lW d 2 3 + 11 e^d^r] 3 + 4e 3 f3 3 d 2 3 {-9 e 2 w 2 d 1 2 d 2 + d 2 3 
+ 13 s 3 d 1 3 w 3 + 11 wdxe d 2 2 )rf + d 2 V/3 4 (113 e 2 w 2 d 2 + d 2 2 + 22 e d lW d 2 )r] - 4 d 2 5 e 6 wd 1 (3 5 ), 

In the following, linear stability analysis yields the bifurcation diagram with r = 0.5, e — 1, (3 — 0.6, B = 0.4, 
r) = 0.25, w = 0.4,d 2 = 1 shown in Fig. [HA). 

The Hopf bifurcation line and the Turing bifurcation curve separate the parametric space into four distinct domains. 
In domain I, located below all two bifurcation lines, the steady state is the only stable solution of the system. Domain 
II is the region of pure Turing instability, while domain III is the region of pure Hopf instability. In domain IV, which 
is located above all two bifurcation lines, both Hopf and Turing instability occur. 

To see the relation between the real and the imaginary parts of the eigenvalue A(fc), we plot in Fig. HJB)~(F) the 
real and the imaginary parts of the eigenvalue at different K with r = 0.5, e — 1, f3 — 0.6, B = 0.4, 77 = 0.25, w = 0.4 
di = 0.01 and d 2 = 1 for the system [5] 

From the definition of Hopf and Turing bifurcation, we know that the relation between the real, the imaginary parts 
of the eigenvalue \{k) determine the bifurcation type. The relation between Re(A(fc)), Im(A(fc)) and k are shown in 
figure QJB)-(F). Figure HJB) illustrate the case of parameter locate in domain I in figure QJA), K = 2.131712170, 
the critical value of Turing bifurcation, in this case, Re(A(fc)) > at k = while Im(A(fc)) 7^ 0. In figure QJC)(D), 
K = 2.2 and K = 2.6, the parameter locate in domain II, the pure Turing instability occurs, one can see that at 
k = 0, Re(A(fc)) = 0, Im(A(») ^ 0. Figure^E), K = 2.631578947, the critical value of Hopf bifurcation, in this case, 
Re(A(fc)) = at k — while Im(A(fe)) 7^ 0. When K = 3.0, parameter locate in domain IV, figure QJF) indicate that 
at k = 0, Re(A(fc)) > 0, Im(A(fc)) ^ 0. 

3. SPATIOTEMPORAL PATTERN FORMATION 

In this section, we perform extensive numerical simulations of the spatially extended model ^ in two-dimensional 
spaces, and the qualitative results are shown here. All our numerical simulations employ the periodic Neumann 
(zero-flux) boundary conditions with a system size of 200 x 200 space units and r = 0.5, e — 1, (3 — 0.6, B — 0.4, 
i] = 0.25, w = 0.4 d% = 0.01 and d 2 = 1. The EqJ2] are solved numerically in two-dimensional space using a finite 
difference approximation for the spatial derivatives and an explicit Euler method for the time integration with a time 
stepsize of At = 0.01 and space stepsize (lattice constant) Ah — 0.25 (see, for details, 0). The scale of the space 
and time are average to the Euler method. The initial density distribution corresponds to random perturbations 
around the stationary state (N* , P*) in model $2$ with a variance significantly lower than the amplitude of the final 
patterns, which seems to be more general from the biological point of view. When the system reached a stable state 
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FIG. 1: (A) K - di Bifurcation diagram for the modeirjwith r = 0.5, e = 1, j3 = 0.6, B = 0.4, 77 = 0.25, w = 0.4 and d 2 = 1. 
Hopf and Turing bifurcation lines separate the parameter space into four domains. And the Hopf- Turing bifurcation point is 
(di,K) — (0.02742,2.63158). In (B) — (F), Re(A(fc)) and Im(A(fc)) are shown by solid curves and dotted curves, respectively. 
The other parameters are: di = 0.01, the bifurcation parameter K: (B) 2.131712170, the critical value of Kt)', (C) 2.2; (D) 
2.6; (E) 2.631578947 the critical value of K H ); (F) 3.0. 

(stationary or oscillatory), we took a snapshot with yellow levels linearly proportional to the free species density and 
red corresponding to high while blue corresponding to low. 

In the numerical simulations, different types of dynamics are observed and we have found that the distributions of 
predator and prey are always of the same type. Consequently, we can restrict our analysis of pattern formation to 
one distribution. In this section, we show the distribution of prey, for instance. 

Fig. [5] shows the evolution of the spatial pattern of prey at 0, 10000, 200000 and 300000 iterations, with random 
small perturbation of the stationary solution (N* , P*) — (0.5094, 0.7829) of the spatially homogeneous systems when 
K = 2.2, located in domain II, slightly more than the Turing bifurcation threshold Kt = 2.1317 and less than the 
Hopf bifurcation threshold Kh = 2.6316. In this case, one can see that for model [2 the random initial distribution 
(c.f., Fig.[^A)) leads to the formation of a regular macroscopic spotted pattern which prevails over the whole domain 
at last, and the dynamics of the system does not undergo any further changes (c.f., Fig.[^D)). 

Fig. [3] shows the evolution of the spatial pattern of prey at 0, 10000, 30000 and 100000 iterations when K = 2.6 
which is more than Kt and slightly less than Kh- Although the dynamics of the system starts from the stationary 
solution (N*,P*) = (0.5252,0.8382), there is an essential difference in the case above. From the snapshots, one can 
see that the steady state of spotted pattern and the stripe-like pattern coexist (c.f., Fig. ElD))- 

In Fig. 21 Kh < K = 3.0, i.e., parameters in domain IV, both Hopf and Turing instability occur. In this case, 
(JV*,P*) = (0.5380,0.8831). One can see that the evolution of the spatial pattern of prey at 0, 20000, 50000 and 
100000 iterations. After the spatial chaos patterns (c.f., Fig. QJB)), a regular stationary stripe- like spatial state 
emerges (c.f., Fig.|4jD)). 

When K = 2.65, d\ = 0.04, i.e., parameters in domain III (c.f., Fig. [HA)), pure Hopf instability occurs. As an 
example, the formation of a regular macroscopic two-dimensional spatial pattern, the spiral pattern, is shown in 
Fig. [5] with a system size of 400 x 400 space units. One can see that for model [21 the random initial distribution 
around the steady state (N* , P*) = (0.5270, 0.8443), a unstable focus of models leads to the formation of the spiral 
pattern in the domain (c.f., Fig.[5fD)). In other words, in this situation, spatially uniform steady-state predator-prey 
coexistence is no longer. Small random fluctuations will be strongly amplified by diffusion, leading to nonuniform 
population distributions. From the analysis in section II, we find with these parameters in domain III, the pattern 
formation, i.e., the spiral pattern, arises from pure Hopf instability. In order to make it clearer, in Fig. [5J we show 
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FIG. 2: Dynamics of the time evolution of the prey of model [2] with di = 0.01, Kt < K = 2.2 < Kh- (A) iteration, (B) 
10000 iterations, (C) 200000 iterations, (D) 300000 iterations. 




FIG. 3: Dynamics of the time evolution of the prey of model [5] with d\ = 0.01, Kt < K = 2.6 < Kh- (A) iteration, (B) 
10000 iterations, (C) 30000 iterations, (D) 100000 iterations. 




FIG. 5: Dynamics of the time evolution of the prey of model [2] with di = 0.04, Kt < Kh < K — 2.65. (A) iteration, (B) 
30000 iterations, (C) 70000 iterations, (D) 100000 iterations. 
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FIG. 6: Dynamical behaviors of model [2] (A) Phase portrait; (B) Time-series plot; (C) Space-time plots corresponding to 
Fig- IHC). The parameters are the same as those in Fig [5] 

phase portrait (Fig. [BJA)), time series plot (Fig. E^B)) and space-time plots (Fig. EJC)), a one-dimensional example 
corresponding to Fig.[5fC). And the method of space-time plots is that let y be a constant (here, y = 200, the center 
line of each snapshots), from each pattern snapshots, choose the line y — 200, and pile these lines in-time-order. The 
space-time plots show the evolution process of the prey N throughout time t and space x. Fig. [6^A) exhibits the 
"local" phase plane of the system obtained in a fixed point (0.5270, 0.8443) inside the region invaded by the irregular 
spatiotemporal oscillations, and the trajectory fills nearly the whole domain inside the limit cycle. Fig.^B) illustrates 
the evolution process of prey density with time, periodic oscillations in time finally. From Fig. [6]JC), one can clearly 
see that a spiral wave emerges, and the system gives rise to uniform oscillations in space and periodic oscillations in 
time. 

Comparing Fig. [2j Fig. [3] with Fig. H Fig. EJ we can see that the bifurcation parameter K determines the type of 
the pattern formation even with the same parameters, e.g., r, e, (3, B, rj, w, d,2- In domain II of Fig. HJA), the closer 
K is to Kt, the more distinct the spotted spatial pattern becomes (c.f., Fig. [2jD)). When K is much closer to Kh, 
the spotted and stripe-like patterns coexist (c.f., Fig.[3][D)). When K is bigger than Kh and far away from the Turing 
bifurcation value Kt, the distinct stripe-like pattern emerges. When K is bigger than Kh and smaller than Kt, the 
spiral wave pattern occurs (c.f., Fig. EJC,D)). So we may draw a conclusion that for model [2] pure Turing instability 
gives birth to the spotted pattern, pure Hopf bifurcation gives birth to the spiral wave pattern, and both of them give 
birth to the stripe-like pattern. 

4. CONCLUSIONS AND REMARKS 

In this paper, we have presented a theoretical analysis of evolutionary processes that involves organisms distribution 
and their interaction of spatially distributed population with local diffusion. And the numerical simulations were 
consistent with the predictions drawn from the bifurcation analysis, i.e., Hopf bifurcation and Turing instability. In 
the domain II of Fig. [TJ A) , the stationary state of periodic spotted pattern exists when the parameters are near the 
Turing bifurcation line, while near the Hopf bifurcation line, both the spotted pattern and the stripe-like pattern 
coexist. When the parameters are located in domain III, pure Hopf instability occurs, and the spiral wave pattern 
emerges. When the parameters are located in domain IV, both Hopf and Turing instability occur, and the stationary 
state of stripe-like pattern exists. 

Turing instability is relevant not only in reaction-diffusion systems, but also in describing other dissipative structures, 
which can be understood in terms of diffusion-driven instability. In addition to the biological relevance of Turing 
systems, their ability to generate structure is of great interest from the point of view of physics. There are various 
physical systems that show similar phenomena, although the underlying mechanisms can be very different. Thus, most 
of the research in the field relies on experiments and numerical simulations justified by an analytical examination. In 
addition, in simulations one may study pattern formation under constraints that are beyond the reach of experiments 
and the numerical data is also easy to analyze. 

In Ref. 2], it's indicated that the basic idea of diffusion-driven instability in a reaction-diffusion model can be 
understood in terms of an activator-inhibitor system. And a random increase of activator species (prey, N) has a 
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positive effect on the creation rate of both activator and inhibitor (predator, P) species. In other words, random 
fluctuations may cause a nonuniform prey density. This elevated prey density has a positive effect both on prey and 
predator population growth rates. Following D. Alonso et al we give the discussion to model([2|). From Eqs.([2]), 
we can obtain the following equations: 

1 dN _ N\ /3P 1 dP _. ef3N 

N dt V K> B+N+wP' P dt ~ B+N+wP '/• V ±u ') 

Similar to Ref. [2J, the first equation in Eqs.(|10p is a one- humped function of prey density, the numerical result can 
be found in [20j | , and prey growth rate can be increased by a higher local prey density at least in a range of parameter 
values. On the other hand, the second equation in Eqs.ljlOJ). i.e., predator numerical response, is an ever- increasing 
function of N, and high prey density always has a positive influence on predator growth. More importantly, inhibitor 
species (predator, P) must diffuse faster than activator species (prey, N), for an increment in inhibitor species may 
have a negative effect on formation rate of both species. Thus, as random fluctuations increase local prey density 
over its equilibrium value, prey population undergoes an accelerated growth. Simultaneously, predator population 
also increases, but as predators diffuse faster than prey, they disperse away from the center of prey outbreaks. If 
relative diffusion (tfe/di) is large enough, prey growth rate will reach negative values and prey population will be 
driven by predators to a very low level in those regions. The final result is the formation of patches of high prey 
density surrounded by areas of low prey density. Predators follow the same pattern. 



In Ref. [ll|], Neuhauser and Pacala formulated the Lotka-Volterra model as a spatial one. They found the striking 



result that the coexistence of patterns is actually harder to get in the spatial model than in the non-spatial one. One 
reason can be traced to how local interaction between individual members of the species are represented in the model. 
In this study, our results show that the predator-prey model with Beddington-DeAngelis-type functional response 
and reaction-diffusion (e.g., Eqs.[2|) also represents rich spatial dynamics, such as spotted pattern, stripe-like pattern, 
coexistence of both stripe-like and spotted pattern, spiral pattern, etc. It will be useful for studying the dynamic 
complexity of ecosystems or physical systems. In Ref. |9|, it is indicated that reaction-diffusion models provide a way 
to translate local assumptions about the movement, mortality, and reproduction of individuals into global conclusions 
about the persistence or extinction of populations and the coexistence of interacting species. They can be derived 
mechanistically via rescaling from models of individual movement which are based on random walks, i.e., small random 
perturbation of the stationary solution (N*,P*) of the spatially homogeneous model ^j. Reaction-diffusion system, 
i.e., model |2|). is spatially explicit and typically incorporate quantities such as dispersal rates, local growth rates, and 
carrying capacities as parameters which may vary with location or time. 

More interesting, when the parameters are located in domain III, pure Hopf instability occurs, and the spiral wave 
pattern emerges (c.f., Fig. [5j Fig. EJC)). To our knowledge, we haven't got any report about one system that has 
spotted, stripe-like and spiral pattern formation meantime. It's well known that spirals and curves are the most 
fascinating clusters to emerge from the predator-prey model. A spiral will form from a wave front when the rabbit 
line (which is leading the front) overlaps the pursuing line of predator. The prey on the extreme end of the line stop 
moving as there are no predator in their immediate vicinity. However the prey and the predator in the center of the 
line continue moving forward. This forms a small trail of prey at one (or both) ends of the front. These prey start 
breeding and the trailing line of prey thickens and attracts the attention of predator at the end of the fox line that 
turn towards this new source of prey. Thus a spiral forms with predator on the inside and prey on the outside. If the 
original overlap of prey occurs at both ends of the line a double spiral will form. Spirals can also form as a prey blob 



collapses after predator eat into it [23|, |24(. Thus, reaction-diffusion system provides a good framework for studying 



questions about the ways that habitat geometry and the size or variation in vital parameters influence population 
dynamics. 

On the other hand, one can see that although there is a small difference between the denominators of the functional 
responses of Michaelis-Menten-type and those of Beddington-DeAngelis-type, there is an enormous gap between them 
in the process of computations. The Turing bifurcation expression of Michaelis-Menten-type predator-prey system is 
simple [iij], while from (J9j> , one can see that the Turing bifurcation analysis requires huge-sized computations, so we 
have to obtain more help via computers. 

In fact, computer aided analysis is useful for nonlinear analysis. And computers have played an important role 
throughout the history of ecology. Today, numerical simulations also play an important role in spatial ecology. There 
are some international mathematical softwares, such as Matlab, Maple, Mathematica, etc. We have finished all our 
symbolic computations in Maple and obtained our pattern snapshots (i.e., numerical simulations) in Matlab for Maple 
is more superior in symbolic computations while Matlab is more superior in numerical computations. 
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